Overall Mortality
Cox regression of entire cohort (not in quartiles)
#Cox regression model
cox_all <- coxph(Surv(survival, censor) ~ Male+Age+non_white+APACHE+Immsupp+Total, data=imgs)
summary(cox_all)
## Call:
## coxph(formula = Surv(survival, censor) ~ Male + Age + non_white +
## APACHE + Immsupp + Total, data = imgs)
##
## n= 560, number of events= 369
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Male1 0.055338 1.056898 0.105766 0.523 0.60082
## Age 0.017286 1.017436 0.003438 5.027 4.97e-07 ***
## non_white1 -0.344413 0.708636 0.141348 -2.437 0.01482 *
## APACHE 0.026503 1.026857 0.006851 3.868 0.00011 ***
## Immsupp1 0.965437 2.625934 0.110705 8.721 < 2e-16 ***
## Total 0.053079 1.054513 0.012027 4.413 1.02e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Male1 1.0569 0.9462 0.8590 1.3004
## Age 1.0174 0.9829 1.0106 1.0243
## non_white1 0.7086 1.4112 0.5372 0.9348
## APACHE 1.0269 0.9738 1.0132 1.0407
## Immsupp1 2.6259 0.3808 2.1137 3.2622
## Total 1.0545 0.9483 1.0299 1.0797
##
## Concordance= 0.695 (se = 0.016 )
## Rsquare= 0.255 (max possible= 1 )
## Likelihood ratio test= 165.1 on 6 df, p=<2e-16
## Wald test = 158.3 on 6 df, p=<2e-16
## Score (logrank) test = 167.7 on 6 df, p=<2e-16
cox_all_phx <- coxph(Surv(survival, censor) ~ Male+Age+non_white+APACHE+Immsupp+Total+copd+chf, data=imgs)
summary(cox_all_phx)
## Call:
## coxph(formula = Surv(survival, censor) ~ Male + Age + non_white +
## APACHE + Immsupp + Total + copd + chf, data = imgs)
##
## n= 559, number of events= 369
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Male1 0.056936 1.058588 0.106170 0.536 0.591771
## Age 0.017627 1.017783 0.003471 5.078 3.81e-07 ***
## non_white1 -0.331374 0.717937 0.141693 -2.339 0.019353 *
## APACHE 0.026032 1.026373 0.006864 3.792 0.000149 ***
## Immsupp1 0.960629 2.613340 0.111177 8.641 < 2e-16 ***
## Total 0.057090 1.058751 0.012197 4.681 2.86e-06 ***
## copd1 0.038081 1.038815 0.145341 0.262 0.793311
## chf1 -0.311076 0.732658 0.201684 -1.542 0.122978
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Male1 1.0586 0.9447 0.8597 1.3035
## Age 1.0178 0.9825 1.0109 1.0247
## non_white1 0.7179 1.3929 0.5438 0.9478
## APACHE 1.0264 0.9743 1.0127 1.0403
## Immsupp1 2.6133 0.3827 2.1017 3.2496
## Total 1.0588 0.9445 1.0337 1.0844
## copd1 1.0388 0.9626 0.7813 1.3812
## chf1 0.7327 1.3649 0.4934 1.0879
##
## Concordance= 0.696 (se = 0.016 )
## Rsquare= 0.259 (max possible= 1 )
## Likelihood ratio test= 167.9 on 8 df, p=<2e-16
## Wald test = 160.4 on 8 df, p=<2e-16
## Score (logrank) test = 169.9 on 8 df, p=<2e-16
cox_ards <- coxph(Surv(survival, censor) ~ Male+Age+non_white+APACHE+Immsupp+Total+copd+chf, data=imgs_ards)
summary(cox_ards)
## Call:
## coxph(formula = Surv(survival, censor) ~ Male + Age + non_white +
## APACHE + Immsupp + Total + copd + chf, data = imgs_ards)
##
## n= 123, number of events= 89
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Male1 0.022694 1.022954 0.231944 0.098 0.9221
## Age 0.016957 1.017102 0.007723 2.196 0.0281 *
## non_white1 -0.408296 0.664782 0.350104 -1.166 0.2435
## APACHE 0.034026 1.034611 0.015078 2.257 0.0240 *
## Immsupp1 1.504820 4.503344 0.249986 6.020 1.75e-09 ***
## Total -0.006605 0.993416 0.032660 -0.202 0.8397
## copd1 0.299967 1.349814 0.354405 0.846 0.3973
## chf1 0.114992 1.121865 0.490740 0.234 0.8147
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Male1 1.0230 0.9776 0.6493 1.612
## Age 1.0171 0.9832 1.0018 1.033
## non_white1 0.6648 1.5043 0.3347 1.320
## APACHE 1.0346 0.9665 1.0045 1.066
## Immsupp1 4.5033 0.2221 2.7590 7.351
## Total 0.9934 1.0066 0.9318 1.059
## copd1 1.3498 0.7408 0.6739 2.704
## chf1 1.1219 0.8914 0.4288 2.935
##
## Concordance= 0.713 (se = 0.033 )
## Rsquare= 0.357 (max possible= 0.998 )
## Likelihood ratio test= 54.41 on 8 df, p=6e-09
## Wald test = 48.82 on 8 df, p=7e-08
## Score (logrank) test = 54.79 on 8 df, p=5e-09
cox_ards_qt <- coxph(Surv(survival, censor) ~ Male+Age+non_white+APACHE+Immsupp+quartile+copd+chf, data=imgs_ards)
summary(cox_ards_qt)
## Call:
## coxph(formula = Surv(survival, censor) ~ Male + Age + non_white +
## APACHE + Immsupp + quartile + copd + chf, data = imgs_ards)
##
## n= 123, number of events= 89
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Male1 0.019116 1.019300 0.233130 0.082 0.9346
## Age 0.017232 1.017382 0.007756 2.222 0.0263 *
## non_white1 -0.407777 0.665127 0.349985 -1.165 0.2440
## APACHE 0.035300 1.035930 0.016508 2.138 0.0325 *
## Immsupp1 1.510938 4.530978 0.253880 5.951 2.66e-09 ***
## quartile2 -0.103852 0.901358 0.325526 -0.319 0.7497
## quartile3 -0.076088 0.926735 0.323918 -0.235 0.8143
## quartile4 -0.119369 0.887481 0.321491 -0.371 0.7104
## copd1 0.294323 1.342217 0.372813 0.789 0.4298
## chf1 0.101239 1.106541 0.493868 0.205 0.8376
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Male1 1.0193 0.9811 0.6454 1.610
## Age 1.0174 0.9829 1.0020 1.033
## non_white1 0.6651 1.5035 0.3350 1.321
## APACHE 1.0359 0.9653 1.0029 1.070
## Immsupp1 4.5310 0.2207 2.7548 7.452
## quartile2 0.9014 1.1094 0.4762 1.706
## quartile3 0.9267 1.0791 0.4912 1.749
## quartile4 0.8875 1.1268 0.4726 1.667
## copd1 1.3422 0.7450 0.6464 2.787
## chf1 1.1065 0.9037 0.4203 2.913
##
## Concordance= 0.713 (se = 0.033 )
## Rsquare= 0.358 (max possible= 0.998 )
## Likelihood ratio test= 54.53 on 10 df, p=4e-08
## Wald test = 48.67 on 10 df, p=5e-07
## Score (logrank) test = 54.85 on 10 df, p=3e-08
schoenfeld_res <- cox.zph(cox_all)
plot(schoenfeld_res)






Cox regression of entire cohort by quartile
#create quartiles
imgs <- imgs %>% mutate(quartile = ntile(Total,4))
imgs$quartile <- as.factor(imgs$quartile)
#Cox regressions
cox_qt <- survfit(Surv(survival, censor) ~ quartile, data=imgs) #for graph
cox_qt
## Call: survfit(formula = Surv(survival, censor) ~ quartile, data = imgs)
##
## n events median 0.95LCL 0.95UCL
## quartile=1 140 64 3117 1802 NA
## quartile=2 140 103 206 136 526
## quartile=3 140 97 296 104 712
## quartile=4 140 105 102 38 244
#cox_quart <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+quartile, data=imgs) #for betas
#summary(cox_quart)
cox_ards_q <- coxph(Surv(survival, censor) ~ Male+Age+non_white+APACHE+Immsupp+quartile, data=imgs_ards)
summary(cox_ards_q)
## Call:
## coxph(formula = Surv(survival, censor) ~ Male + Age + non_white +
## APACHE + Immsupp + quartile, data = imgs_ards)
##
## n= 123, number of events= 89
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Male1 -0.006831 0.993193 0.221369 -0.031 0.9754
## Age 0.017974 1.018136 0.007631 2.356 0.0185 *
## non_white1 -0.410447 0.663354 0.347343 -1.182 0.2373
## APACHE 0.035740 1.036387 0.016526 2.163 0.0306 *
## Immsupp1 1.478997 4.388542 0.249134 5.937 2.91e-09 ***
## quartile2 -0.119471 0.887390 0.323193 -0.370 0.7116
## quartile3 -0.010618 0.989438 0.309665 -0.034 0.9726
## quartile4 -0.111176 0.894781 0.318036 -0.350 0.7267
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Male1 0.9932 1.0069 0.6436 1.533
## Age 1.0181 0.9822 1.0030 1.033
## non_white1 0.6634 1.5075 0.3358 1.310
## APACHE 1.0364 0.9649 1.0034 1.071
## Immsupp1 4.3885 0.2279 2.6931 7.151
## quartile2 0.8874 1.1269 0.4710 1.672
## quartile3 0.9894 1.0107 0.5393 1.815
## quartile4 0.8948 1.1176 0.4797 1.669
##
## Concordance= 0.712 (se = 0.033 )
## Rsquare= 0.355 (max possible= 0.998 )
## Likelihood ratio test= 53.91 on 8 df, p=7e-09
## Wald test = 48.38 on 8 df, p=8e-08
## Score (logrank) test = 54.41 on 8 df, p=6e-09
cox_ards_qp <- survfit(Surv(survival, censor) ~ quartile, data=imgs_ards)
plot_ards <- ggsurvplot(cox_ards_qp, data=imgs_ards, pval=TRUE, pval.size=3,
ggtheme = theme_minimal(), legend.title="Quartile", legend.labs=c("Q1","Q2","Q3","Q4"),
risk.table=TRUE, tables.height = 0.2, tables.theme = theme_cleantable(),
xlim=c(0,3650), break.time.by=365.25, xscale="d_m", ylab="",
xlab="Months", font.legend=8, title="Entire cohort", font.title=10,
legend=c(1.1,0.6), risk.table.y.text=FALSE,
risk.table.fontsize=3, risk.table.title="")
plot_ards$plot <- plot_ards$plot + theme(plot.title = element_text(hjust = 0.5))
plot_ards

#Adding copd/chf did not explode the std errors.
contrasts(imgs$quartile) <- contr.treatment(4, base=4)
cox_qt_base4 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs)
summary(cox_qt_base4)
## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp +
## non_white + APACHE + copd + chf + quartile, data = imgs)
##
## n= 559, number of events= 369
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Age 0.015858 1.015984 0.003519 4.506 6.60e-06 ***
## Male1 0.054028 1.055514 0.106463 0.507 0.611817
## Immsupp1 0.995257 2.705419 0.111847 8.898 < 2e-16 ***
## non_white1 -0.347144 0.706703 0.141602 -2.452 0.014224 *
## APACHE 0.025299 1.025622 0.006914 3.659 0.000253 ***
## copd1 0.052707 1.054120 0.147045 0.358 0.720015
## chf1 -0.261249 0.770089 0.201493 -1.297 0.194780
## quartile1 -0.781829 0.457568 0.163799 -4.773 1.81e-06 ***
## quartile2 -0.059434 0.942298 0.142399 -0.417 0.676403
## quartile3 -0.105750 0.899650 0.143789 -0.735 0.462063
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Age 1.0160 0.9843 1.0090 1.0230
## Male1 1.0555 0.9474 0.8567 1.3004
## Immsupp1 2.7054 0.3696 2.1729 3.3685
## non_white1 0.7067 1.4150 0.5354 0.9328
## APACHE 1.0256 0.9750 1.0118 1.0396
## copd1 1.0541 0.9487 0.7902 1.4062
## chf1 0.7701 1.2986 0.5188 1.1430
## quartile1 0.4576 2.1855 0.3319 0.6308
## quartile2 0.9423 1.0612 0.7128 1.2457
## quartile3 0.8996 1.1115 0.6787 1.1925
##
## Concordance= 0.701 (se = 0.016 )
## Rsquare= 0.271 (max possible= 1 )
## Likelihood ratio test= 176.4 on 10 df, p=<2e-16
## Wald test = 162.9 on 10 df, p=<2e-16
## Score (logrank) test = 173.4 on 10 df, p=<2e-16
#Result: Only Q1 is significant
contrasts(imgs$quartile) <- contr.treatment(4, base=3)
cox_qt_base3 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs)
summary(cox_qt_base3)
## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp +
## non_white + APACHE + copd + chf + quartile, data = imgs)
##
## n= 559, number of events= 369
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Age 0.015858 1.015984 0.003519 4.506 6.60e-06 ***
## Male1 0.054028 1.055514 0.106463 0.507 0.611817
## Immsupp1 0.995257 2.705419 0.111847 8.898 < 2e-16 ***
## non_white1 -0.347144 0.706703 0.141602 -2.452 0.014224 *
## APACHE 0.025299 1.025622 0.006914 3.659 0.000253 ***
## copd1 0.052707 1.054120 0.147045 0.358 0.720015
## chf1 -0.261249 0.770089 0.201493 -1.297 0.194780
## quartile1 -0.676079 0.508607 0.165754 -4.079 4.53e-05 ***
## quartile2 0.046316 1.047405 0.147097 0.315 0.752863
## quartile4 0.105750 1.111544 0.143789 0.735 0.462063
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Age 1.0160 0.9843 1.0090 1.0230
## Male1 1.0555 0.9474 0.8567 1.3004
## Immsupp1 2.7054 0.3696 2.1729 3.3685
## non_white1 0.7067 1.4150 0.5354 0.9328
## APACHE 1.0256 0.9750 1.0118 1.0396
## copd1 1.0541 0.9487 0.7902 1.4062
## chf1 0.7701 1.2986 0.5188 1.1430
## quartile1 0.5086 1.9662 0.3675 0.7038
## quartile2 1.0474 0.9547 0.7851 1.3974
## quartile4 1.1115 0.8996 0.8386 1.4734
##
## Concordance= 0.701 (se = 0.016 )
## Rsquare= 0.271 (max possible= 1 )
## Likelihood ratio test= 176.4 on 10 df, p=<2e-16
## Wald test = 162.9 on 10 df, p=<2e-16
## Score (logrank) test = 173.4 on 10 df, p=<2e-16
#Result: Only Q1 is significant
contrasts(imgs$quartile) <- contr.treatment(4, base=2)
cox_qt_base2 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs)
summary(cox_qt_base2)
## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp +
## non_white + APACHE + copd + chf + quartile, data = imgs)
##
## n= 559, number of events= 369
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Age 0.015858 1.015984 0.003519 4.506 6.60e-06 ***
## Male1 0.054028 1.055514 0.106463 0.507 0.611817
## Immsupp1 0.995257 2.705419 0.111847 8.898 < 2e-16 ***
## non_white1 -0.347144 0.706703 0.141602 -2.452 0.014224 *
## APACHE 0.025299 1.025622 0.006914 3.659 0.000253 ***
## copd1 0.052707 1.054120 0.147045 0.358 0.720015
## chf1 -0.261249 0.770089 0.201493 -1.297 0.194780
## quartile1 -0.722395 0.485588 0.163848 -4.409 1.04e-05 ***
## quartile3 -0.046316 0.954740 0.147097 -0.315 0.752863
## quartile4 0.059434 1.061236 0.142399 0.417 0.676403
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Age 1.0160 0.9843 1.0090 1.0230
## Male1 1.0555 0.9474 0.8567 1.3004
## Immsupp1 2.7054 0.3696 2.1729 3.3685
## non_white1 0.7067 1.4150 0.5354 0.9328
## APACHE 1.0256 0.9750 1.0118 1.0396
## copd1 1.0541 0.9487 0.7902 1.4062
## chf1 0.7701 1.2986 0.5188 1.1430
## quartile1 0.4856 2.0594 0.3522 0.6695
## quartile3 0.9547 1.0474 0.7156 1.2738
## quartile4 1.0612 0.9423 0.8028 1.4029
##
## Concordance= 0.701 (se = 0.016 )
## Rsquare= 0.271 (max possible= 1 )
## Likelihood ratio test= 176.4 on 10 df, p=<2e-16
## Wald test = 162.9 on 10 df, p=<2e-16
## Score (logrank) test = 173.4 on 10 df, p=<2e-16
#Result: Only Q1 is significant
#reset factor to base 1
contrasts(imgs$quartile) <- contr.treatment(4, base=1)
cox_qt_base1 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs)
summary(cox_qt_base1)
## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp +
## non_white + APACHE + copd + chf + quartile, data = imgs)
##
## n= 559, number of events= 369
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Age 0.015858 1.015984 0.003519 4.506 6.60e-06 ***
## Male1 0.054028 1.055514 0.106463 0.507 0.611817
## Immsupp1 0.995257 2.705419 0.111847 8.898 < 2e-16 ***
## non_white1 -0.347144 0.706703 0.141602 -2.452 0.014224 *
## APACHE 0.025299 1.025622 0.006914 3.659 0.000253 ***
## copd1 0.052707 1.054120 0.147045 0.358 0.720015
## chf1 -0.261249 0.770089 0.201493 -1.297 0.194780
## quartile2 0.722395 2.059360 0.163848 4.409 1.04e-05 ***
## quartile3 0.676079 1.966154 0.165754 4.079 4.53e-05 ***
## quartile4 0.781829 2.185466 0.163799 4.773 1.81e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Age 1.0160 0.9843 1.0090 1.0230
## Male1 1.0555 0.9474 0.8567 1.3004
## Immsupp1 2.7054 0.3696 2.1729 3.3685
## non_white1 0.7067 1.4150 0.5354 0.9328
## APACHE 1.0256 0.9750 1.0118 1.0396
## copd1 1.0541 0.9487 0.7902 1.4062
## chf1 0.7701 1.2986 0.5188 1.1430
## quartile2 2.0594 0.4856 1.4937 2.8392
## quartile3 1.9662 0.5086 1.4208 2.7209
## quartile4 2.1855 0.4576 1.5853 3.0128
##
## Concordance= 0.701 (se = 0.016 )
## Rsquare= 0.271 (max possible= 1 )
## Likelihood ratio test= 176.4 on 10 df, p=<2e-16
## Wald test = 162.9 on 10 df, p=<2e-16
## Score (logrank) test = 173.4 on 10 df, p=<2e-16
plot_all <- ggsurvplot(cox_qt, data=imgs, pval=TRUE, pval.size=3,
ggtheme = theme_minimal(), legend.title="Quartile", legend.labs=c("Q1","Q2","Q3","Q4"),
risk.table=TRUE, tables.height = 0.2, tables.theme = theme_cleantable(),
xlim=c(0,3650), break.time.by=365.25, xscale="d_m", ylab="",
xlab="Months", font.legend=8, title="Entire cohort", font.title=10,
legend=c(1.1,0.6), risk.table.y.text=FALSE,
risk.table.fontsize=3, risk.table.title="")
plot_all$plot <- plot_all$plot + theme(plot.title = element_text(hjust = 0.5))
plot_all

Cox regression for non-ARDS subset
#non-ARDS patients
imgs_non <- imgs %>% filter(Final_Dx != "ARDS" & Final_Dx != "Sepsis/ARDS") %>% mutate(quartile = as.factor(ntile(Total,4)))
cox_all_non <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+Total, data=imgs_non) #for betas
summary(cox_all_non)
## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp +
## non_white + APACHE + copd + chf + Total, data = imgs_non)
##
## n= 436, number of events= 280
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Age 0.019080 1.019263 0.003961 4.816 1.46e-06 ***
## Male1 0.082921 1.086456 0.123718 0.670 0.50270
## Immsupp1 0.744200 2.104758 0.129506 5.746 9.11e-09 ***
## non_white1 -0.409955 0.663680 0.157584 -2.602 0.00928 **
## APACHE 0.017896 1.018057 0.008259 2.167 0.03025 *
## copd1 0.030619 1.031093 0.160688 0.191 0.84888
## chf1 -0.429238 0.651005 0.227067 -1.890 0.05871 .
## Total 0.068074 1.070444 0.014178 4.801 1.58e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Age 1.0193 0.9811 1.0114 1.0272
## Male1 1.0865 0.9204 0.8525 1.3846
## Immsupp1 2.1048 0.4751 1.6329 2.7129
## non_white1 0.6637 1.5067 0.4873 0.9038
## APACHE 1.0181 0.9823 1.0017 1.0347
## copd1 1.0311 0.9698 0.7525 1.4128
## chf1 0.6510 1.5361 0.4172 1.0159
## Total 1.0704 0.9342 1.0411 1.1006
##
## Concordance= 0.68 (se = 0.018 )
## Rsquare= 0.225 (max possible= 0.999 )
## Likelihood ratio test= 111.3 on 8 df, p=<2e-16
## Wald test = 106.5 on 8 df, p=<2e-16
## Score (logrank) test = 111.1 on 8 df, p=<2e-16
#std errors checked -> ok
cox_all_non_strata <- survfit(Surv(survival, censor) ~ quartile, data = imgs_non) #for KM
#std errors checked ->okay
contrasts(imgs_non$quartile) <- contr.treatment(4, base=4)
cox_non_base4 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs_non)
summary(cox_non_base4)
## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp +
## non_white + APACHE + copd + chf + quartile, data = imgs_non)
##
## n= 436, number of events= 280
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Age 0.017938 1.018100 0.003939 4.553 5.28e-06 ***
## Male1 0.083598 1.087192 0.123690 0.676 0.49912
## Immsupp1 0.791019 2.205644 0.130378 6.067 1.30e-09 ***
## non_white1 -0.416692 0.659224 0.157258 -2.650 0.00806 **
## APACHE 0.015815 1.015940 0.008309 1.903 0.05698 .
## copd1 0.023779 1.024064 0.163289 0.146 0.88422
## chf1 -0.396363 0.672763 0.229551 -1.727 0.08422 .
## quartile1 -0.918378 0.399166 0.191878 -4.786 1.70e-06 ***
## quartile2 -0.330150 0.718816 0.170259 -1.939 0.05249 .
## quartile3 -0.079902 0.923207 0.161820 -0.494 0.62147
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Age 1.0181 0.9822 1.0103 1.0260
## Male1 1.0872 0.9198 0.8531 1.3855
## Immsupp1 2.2056 0.4534 1.7083 2.8478
## non_white1 0.6592 1.5169 0.4844 0.8972
## APACHE 1.0159 0.9843 0.9995 1.0326
## copd1 1.0241 0.9765 0.7436 1.4103
## chf1 0.6728 1.4864 0.4290 1.0550
## quartile1 0.3992 2.5052 0.2740 0.5814
## quartile2 0.7188 1.3912 0.5149 1.0036
## quartile3 0.9232 1.0832 0.6723 1.2678
##
## Concordance= 0.684 (se = 0.018 )
## Rsquare= 0.238 (max possible= 0.999 )
## Likelihood ratio test= 118.4 on 10 df, p=<2e-16
## Wald test = 109.5 on 10 df, p=<2e-16
## Score (logrank) test = 115.3 on 10 df, p=<2e-16
#Result: Only Q1 is significant
contrasts(imgs_non$quartile) <- contr.treatment(4, base=3)
cox_non_base3 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs_non)
summary(cox_non_base3)
## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp +
## non_white + APACHE + copd + chf + quartile, data = imgs_non)
##
## n= 436, number of events= 280
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Age 0.017938 1.018100 0.003939 4.553 5.28e-06 ***
## Male1 0.083598 1.087192 0.123690 0.676 0.49912
## Immsupp1 0.791019 2.205644 0.130378 6.067 1.30e-09 ***
## non_white1 -0.416692 0.659224 0.157258 -2.650 0.00806 **
## APACHE 0.015815 1.015940 0.008309 1.903 0.05698 .
## copd1 0.023779 1.024064 0.163289 0.146 0.88422
## chf1 -0.396363 0.672763 0.229551 -1.727 0.08422 .
## quartile1 -0.838476 0.432369 0.186605 -4.493 7.01e-06 ***
## quartile2 -0.250249 0.778607 0.166932 -1.499 0.13385
## quartile4 0.079902 1.083180 0.161820 0.494 0.62147
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Age 1.0181 0.9822 1.0103 1.0260
## Male1 1.0872 0.9198 0.8531 1.3855
## Immsupp1 2.2056 0.4534 1.7083 2.8478
## non_white1 0.6592 1.5169 0.4844 0.8972
## APACHE 1.0159 0.9843 0.9995 1.0326
## copd1 1.0241 0.9765 0.7436 1.4103
## chf1 0.6728 1.4864 0.4290 1.0550
## quartile1 0.4324 2.3128 0.2999 0.6233
## quartile2 0.7786 1.2843 0.5613 1.0800
## quartile4 1.0832 0.9232 0.7888 1.4874
##
## Concordance= 0.684 (se = 0.018 )
## Rsquare= 0.238 (max possible= 0.999 )
## Likelihood ratio test= 118.4 on 10 df, p=<2e-16
## Wald test = 109.5 on 10 df, p=<2e-16
## Score (logrank) test = 115.3 on 10 df, p=<2e-16
#Result: Only Q1 is significant
contrasts(imgs_non$quartile) <- contr.treatment(4, base=2)
cox_non_base2 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs_non)
summary(cox_non_base2)
## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp +
## non_white + APACHE + copd + chf + quartile, data = imgs_non)
##
## n= 436, number of events= 280
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Age 0.017938 1.018100 0.003939 4.553 5.28e-06 ***
## Male1 0.083598 1.087192 0.123690 0.676 0.49912
## Immsupp1 0.791019 2.205644 0.130378 6.067 1.30e-09 ***
## non_white1 -0.416692 0.659224 0.157258 -2.650 0.00806 **
## APACHE 0.015815 1.015940 0.008309 1.903 0.05698 .
## copd1 0.023779 1.024064 0.163289 0.146 0.88422
## chf1 -0.396363 0.672763 0.229551 -1.727 0.08422 .
## quartile1 -0.588228 0.555311 0.194310 -3.027 0.00247 **
## quartile3 0.250249 1.284345 0.166932 1.499 0.13385
## quartile4 0.330150 1.391177 0.170259 1.939 0.05249 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Age 1.0181 0.9822 1.0103 1.0260
## Male1 1.0872 0.9198 0.8531 1.3855
## Immsupp1 2.2056 0.4534 1.7083 2.8478
## non_white1 0.6592 1.5169 0.4844 0.8972
## APACHE 1.0159 0.9843 0.9995 1.0326
## copd1 1.0241 0.9765 0.7436 1.4103
## chf1 0.6728 1.4864 0.4290 1.0550
## quartile1 0.5553 1.8008 0.3794 0.8127
## quartile3 1.2843 0.7786 0.9260 1.7815
## quartile4 1.3912 0.7188 0.9965 1.9423
##
## Concordance= 0.684 (se = 0.018 )
## Rsquare= 0.238 (max possible= 0.999 )
## Likelihood ratio test= 118.4 on 10 df, p=<2e-16
## Wald test = 109.5 on 10 df, p=<2e-16
## Score (logrank) test = 115.3 on 10 df, p=<2e-16
#Result: Only Q1 is significant
#reset factor to base 1
contrasts(imgs_non$quartile) <- contr.treatment(4, base=1)
cox_non_base1 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs_non)
summary(cox_non_base1)
## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp +
## non_white + APACHE + copd + chf + quartile, data = imgs_non)
##
## n= 436, number of events= 280
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Age 0.017938 1.018100 0.003939 4.553 5.28e-06 ***
## Male1 0.083598 1.087192 0.123690 0.676 0.49912
## Immsupp1 0.791019 2.205644 0.130378 6.067 1.30e-09 ***
## non_white1 -0.416692 0.659224 0.157258 -2.650 0.00806 **
## APACHE 0.015815 1.015940 0.008309 1.903 0.05698 .
## copd1 0.023779 1.024064 0.163289 0.146 0.88422
## chf1 -0.396363 0.672763 0.229551 -1.727 0.08422 .
## quartile2 0.588228 1.800794 0.194310 3.027 0.00247 **
## quartile3 0.838476 2.312840 0.186605 4.493 7.01e-06 ***
## quartile4 0.918378 2.505223 0.191878 4.786 1.70e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Age 1.0181 0.9822 1.0103 1.0260
## Male1 1.0872 0.9198 0.8531 1.3855
## Immsupp1 2.2056 0.4534 1.7083 2.8478
## non_white1 0.6592 1.5169 0.4844 0.8972
## APACHE 1.0159 0.9843 0.9995 1.0326
## copd1 1.0241 0.9765 0.7436 1.4103
## chf1 0.6728 1.4864 0.4290 1.0550
## quartile2 1.8008 0.5553 1.2305 2.6355
## quartile3 2.3128 0.4324 1.6044 3.3341
## quartile4 2.5052 0.3992 1.7200 3.6490
##
## Concordance= 0.684 (se = 0.018 )
## Rsquare= 0.238 (max possible= 0.999 )
## Likelihood ratio test= 118.4 on 10 df, p=<2e-16
## Wald test = 109.5 on 10 df, p=<2e-16
## Score (logrank) test = 115.3 on 10 df, p=<2e-16
plot_non <- ggsurvplot(cox_all_non_strata, data=imgs_non, pval=TRUE, pval.size=3,
ggtheme = theme_minimal(), legend.title="Quartile",
risk.table=TRUE, tables.height = 0.2, tables.theme = theme_cleantable(),
xlim=c(0,3650), break.time.by=365.25, xscale="d_m", ylab="",
xlab="Months", font.legend=8, title="Non-ARDS patients", font.title=10,
legend="none", risk.table.y.text=FALSE,
risk.table.fontsize=3, risk.table.title="")
plot_non$plot <- plot_non$plot + theme(axis.text.y = element_blank(), plot.title = element_text(hjust = 0.5))
plot_non

KM curves side by side
KM_curves <- list()
KM_curves[[1]] <- plot_all
KM_curves[[2]] <- plot_non
arrange_ggsurvplots(KM_curves, print=TRUE, ncol=2, nrow=1)

#std errors checked ->okay
hosp_mort <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs)
summary(hosp_mort)
##
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp +
## APACHE + non_white + copd + chf, family = binomial(link = "logit"),
## data = imgs)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5878 -0.7771 -0.4981 0.8897 2.5633
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.308019 0.619638 -6.952 3.59e-12 ***
## quartile2 0.977419 0.350667 2.787 0.005315 **
## quartile3 1.273048 0.344812 3.692 0.000222 ***
## quartile4 1.313196 0.342296 3.836 0.000125 ***
## Male1 -0.007186 0.211581 -0.034 0.972908
## Age 0.007612 0.006892 1.104 0.269390
## Immsupp1 1.284147 0.215892 5.948 2.71e-09 ***
## APACHE 0.058885 0.013126 4.486 7.25e-06 ***
## non_white1 -0.533414 0.294295 -1.813 0.069907 .
## copd1 -0.132786 0.302055 -0.440 0.660220
## chf1 -0.728473 0.417400 -1.745 0.080939 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 667.55 on 558 degrees of freedom
## Residual deviance: 556.32 on 548 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 578.32
##
## Number of Fisher Scoring iterations: 5
cbind(OR = exp(coef(hosp_mort)),
exp(confint(hosp_mort)),
p = coef(summary(hosp_mort))[,4])
## OR 2.5 % 97.5 % p
## (Intercept) 0.01346018 0.003827569 0.04361803 3.589179e-12
## quartile2 2.65758874 1.357107104 5.40098142 5.314693e-03
## quartile3 3.57172343 1.848338033 7.18684749 2.224931e-04
## quartile4 3.71803703 1.934507631 7.44744932 1.248334e-04
## Male1 0.99284009 0.655723601 1.50471217 9.729076e-01
## Age 1.00764099 0.994182653 1.02146055 2.693897e-01
## Immsupp1 3.61158417 2.374076826 5.54060270 2.712762e-09
## APACHE 1.06065325 1.034088295 1.08880700 7.247657e-06
## non_white1 0.58659886 0.322492002 1.02783849 6.990708e-02
## copd1 0.87565251 0.477472759 1.56698187 6.602203e-01
## chf1 0.48264533 0.202196883 1.05496607 8.093892e-02
contrasts(imgs$quartile) <- contr.treatment(4, base=2)
inhosp_quart_base2 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs)
summary(inhosp_quart_base2)
##
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp +
## APACHE + non_white + copd + chf, family = binomial(link = "logit"),
## data = imgs)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5878 -0.7771 -0.4981 0.8897 2.5633
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.330600 0.618051 -5.389 7.09e-08 ***
## quartile1 -0.977419 0.350667 -2.787 0.00531 **
## quartile3 0.295629 0.283931 1.041 0.29778
## quartile4 0.335777 0.279645 1.201 0.22986
## Male1 -0.007186 0.211581 -0.034 0.97291
## Age 0.007612 0.006892 1.104 0.26939
## Immsupp1 1.284147 0.215892 5.948 2.71e-09 ***
## APACHE 0.058885 0.013126 4.486 7.25e-06 ***
## non_white1 -0.533414 0.294295 -1.813 0.06991 .
## copd1 -0.132786 0.302055 -0.440 0.66022
## chf1 -0.728473 0.417400 -1.745 0.08094 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 667.55 on 558 degrees of freedom
## Residual deviance: 556.32 on 548 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 578.32
##
## Number of Fisher Scoring iterations: 5
#result: Q1is significant, 3 and 4 are not
cbind(OR = exp(coef(inhosp_quart_base2)),
exp(confint(inhosp_quart_base2)),
p = coef(summary(inhosp_quart_base2))[,4])
## OR 2.5 % 97.5 % p
## (Intercept) 0.03577163 0.01031823 0.1168662 7.089952e-08
## quartile1 0.37628094 0.18515153 0.7368615 5.314693e-03
## quartile3 1.34397147 0.77105057 2.3519011 2.977824e-01
## quartile4 1.39902649 0.80976693 2.4284928 2.298585e-01
## Male1 0.99284009 0.65572360 1.5047122 9.729076e-01
## Age 1.00764099 0.99418265 1.0214605 2.693897e-01
## Immsupp1 3.61158417 2.37407683 5.5406027 2.712762e-09
## APACHE 1.06065325 1.03408829 1.0888070 7.247657e-06
## non_white1 0.58659886 0.32249200 1.0278385 6.990708e-02
## copd1 0.87565251 0.47747276 1.5669819 6.602203e-01
## chf1 0.48264533 0.20219688 1.0549661 8.093892e-02
contrasts(imgs$quartile) <- contr.treatment(4, base=3)
inhosp_quart_base3 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs)
summary(inhosp_quart_base3)
##
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp +
## APACHE + non_white + copd + chf, family = binomial(link = "logit"),
## data = imgs)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5878 -0.7771 -0.4981 0.8897 2.5633
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.034971 0.599423 -5.063 4.12e-07 ***
## quartile1 -1.273048 0.344812 -3.692 0.000222 ***
## quartile2 -0.295629 0.283931 -1.041 0.297782
## quartile4 0.040148 0.271342 0.148 0.882375
## Male1 -0.007186 0.211581 -0.034 0.972908
## Age 0.007612 0.006892 1.104 0.269390
## Immsupp1 1.284147 0.215892 5.948 2.71e-09 ***
## APACHE 0.058885 0.013126 4.486 7.25e-06 ***
## non_white1 -0.533414 0.294295 -1.813 0.069907 .
## copd1 -0.132786 0.302055 -0.440 0.660220
## chf1 -0.728473 0.417400 -1.745 0.080939 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 667.55 on 558 degrees of freedom
## Residual deviance: 556.32 on 548 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 578.32
##
## Number of Fisher Scoring iterations: 5
#result: Q1 is significant, 2 and 4 are not
cbind(OR = exp(coef(inhosp_quart_base3)),
exp(confint(inhosp_quart_base3)),
p = coef(summary(inhosp_quart_base3))[,4])
## OR 2.5 % 97.5 % p
## (Intercept) 0.04807605 0.01441297 0.1517703 4.123787e-07
## quartile1 0.27997689 0.13914307 0.5410266 2.224931e-04
## quartile2 0.74406342 0.42518794 1.2969318 2.977824e-01
## quartile4 1.04096443 0.61117357 1.7737158 8.823749e-01
## Male1 0.99284009 0.65572360 1.5047122 9.729076e-01
## Age 1.00764099 0.99418265 1.0214605 2.693897e-01
## Immsupp1 3.61158417 2.37407683 5.5406027 2.712762e-09
## APACHE 1.06065325 1.03408829 1.0888070 7.247657e-06
## non_white1 0.58659886 0.32249200 1.0278385 6.990708e-02
## copd1 0.87565251 0.47747276 1.5669819 6.602203e-01
## chf1 0.48264533 0.20219688 1.0549661 8.093892e-02
contrasts(imgs$quartile) <- contr.treatment(4, base=4)
inhosp_quart_base4 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs)
summary(inhosp_quart_base4)
##
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp +
## APACHE + non_white + copd + chf, family = binomial(link = "logit"),
## data = imgs)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5878 -0.7771 -0.4981 0.8897 2.5633
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.994823 0.608013 -4.926 8.41e-07 ***
## quartile1 -1.313196 0.342296 -3.836 0.000125 ***
## quartile2 -0.335777 0.279645 -1.201 0.229858
## quartile3 -0.040148 0.271342 -0.148 0.882375
## Male1 -0.007186 0.211581 -0.034 0.972908
## Age 0.007612 0.006892 1.104 0.269390
## Immsupp1 1.284147 0.215892 5.948 2.71e-09 ***
## APACHE 0.058885 0.013126 4.486 7.25e-06 ***
## non_white1 -0.533414 0.294295 -1.813 0.069907 .
## copd1 -0.132786 0.302055 -0.440 0.660220
## chf1 -0.728473 0.417400 -1.745 0.080939 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 667.55 on 558 degrees of freedom
## Residual deviance: 556.32 on 548 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 578.32
##
## Number of Fisher Scoring iterations: 5
#result: Q1 is significant, 2 and 3 are not
cbind(OR = exp(coef(inhosp_quart_base4)),
exp(confint(inhosp_quart_base4)),
p = coef(summary(inhosp_quart_base4))[,4])
## OR 2.5 % 97.5 % p
## (Intercept) 0.05004546 0.01474234 0.1605245 8.410630e-07
## quartile1 0.26895913 0.13427416 0.5169274 1.248334e-04
## quartile2 0.71478275 0.41177804 1.2349232 2.298585e-01
## quartile3 0.96064762 0.56378818 1.6361964 8.823749e-01
## Male1 0.99284009 0.65572360 1.5047122 9.729076e-01
## Age 1.00764099 0.99418265 1.0214605 2.693897e-01
## Immsupp1 3.61158417 2.37407683 5.5406027 2.712762e-09
## APACHE 1.06065325 1.03408829 1.0888070 7.247657e-06
## non_white1 0.58659886 0.32249200 1.0278385 6.990708e-02
## copd1 0.87565251 0.47747276 1.5669819 6.602203e-01
## chf1 0.48264533 0.20219688 1.0549661 8.093892e-02
#reset factor to base 1
contrasts(imgs$quartile) <- contr.treatment(4, base=1)
#std errors checked ->okay
hosp_mort_non_r1 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs_non)
summary(hosp_mort_non_r1)
##
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp +
## APACHE + non_white + copd + chf, family = binomial(link = "logit"),
## data = imgs_non)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4147 -0.6999 -0.4993 -0.2650 2.5139
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.656926 0.695942 -5.255 1.48e-07 ***
## quartile2 0.422174 0.414174 1.019 0.308052
## quartile3 1.041093 0.382233 2.724 0.006455 **
## quartile4 0.863732 0.393273 2.196 0.028073 *
## Male1 -0.114983 0.253055 -0.454 0.649556
## Age 0.007204 0.008351 0.863 0.388323
## Immsupp1 0.918398 0.259996 3.532 0.000412 ***
## APACHE 0.047948 0.016177 2.964 0.003037 **
## non_white1 -0.843357 0.373297 -2.259 0.023870 *
## copd1 -0.123341 0.346305 -0.356 0.721718
## chf1 -0.862727 0.525182 -1.643 0.100441
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 451.96 on 435 degrees of freedom
## Residual deviance: 398.98 on 425 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 420.98
##
## Number of Fisher Scoring iterations: 5
cbind(OR = exp(coef(hosp_mort_non_r1)),
exp(confint(hosp_mort_non_r1)),
p = coef(summary(hosp_mort_non_r1))[,4])
## OR 2.5 % 97.5 % p
## (Intercept) 0.02581175 0.006256306 0.09638051 1.483154e-07
## quartile2 1.52527456 0.682180953 3.49756666 3.080525e-01
## quartile3 2.83231015 1.363691206 6.15708813 6.455294e-03
## quartile4 2.37199718 1.113636707 5.25176753 2.807291e-02
## Male1 0.89138126 0.542383131 1.46589497 6.495556e-01
## Age 1.00723006 0.991001173 1.02407737 3.883234e-01
## Immsupp1 2.50527300 1.507060048 4.18531889 4.118774e-04
## APACHE 1.04911626 1.016565867 1.08333341 3.036661e-03
## non_white1 0.43026373 0.196597525 0.86158480 2.387012e-02
## copd1 0.88396182 0.437157094 1.71182926 7.217176e-01
## chf1 0.42200962 0.134620177 1.09600071 1.004410e-01
#ref level =2
contrasts(imgs_non$quartile) <- contr.treatment(4, base = 2)
hosp_mort_non_r2 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs_non)
summary(hosp_mort_non_r2)
##
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp +
## APACHE + non_white + copd + chf, family = binomial(link = "logit"),
## data = imgs_non)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4147 -0.6999 -0.4993 -0.2650 2.5139
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.234751 0.738440 -4.381 1.18e-05 ***
## quartile1 -0.422174 0.414174 -1.019 0.308052
## quartile3 0.618918 0.349294 1.772 0.076409 .
## quartile4 0.441558 0.358730 1.231 0.218363
## Male1 -0.114983 0.253055 -0.454 0.649556
## Age 0.007204 0.008351 0.863 0.388323
## Immsupp1 0.918398 0.259996 3.532 0.000412 ***
## APACHE 0.047948 0.016177 2.964 0.003037 **
## non_white1 -0.843357 0.373297 -2.259 0.023870 *
## copd1 -0.123341 0.346305 -0.356 0.721718
## chf1 -0.862727 0.525182 -1.643 0.100441
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 451.96 on 435 degrees of freedom
## Residual deviance: 398.98 on 425 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 420.98
##
## Number of Fisher Scoring iterations: 5
cbind(OR = exp(coef(hosp_mort_non_r2)),
exp(confint(hosp_mort_non_r2)),
p = coef(summary(hosp_mort_non_r2))[,4])
## OR 2.5 % 97.5 % p
## (Intercept) 0.0393700 0.008820534 0.1606021 1.183977e-05
## quartile1 0.6556197 0.285913064 1.4658867 3.080525e-01
## quartile3 1.8569182 0.942948850 3.7287019 7.640892e-02
## quartile4 1.5551280 0.773205563 3.1735108 2.183629e-01
## Male1 0.8913813 0.542383131 1.4658950 6.495556e-01
## Age 1.0072301 0.991001173 1.0240774 3.883234e-01
## Immsupp1 2.5052730 1.507060048 4.1853189 4.118774e-04
## APACHE 1.0491163 1.016565867 1.0833334 3.036661e-03
## non_white1 0.4302637 0.196597525 0.8615848 2.387012e-02
## copd1 0.8839618 0.437157094 1.7118293 7.217176e-01
## chf1 0.4220096 0.134620177 1.0960007 1.004410e-01
#ref level = 3
contrasts(imgs_non$quartile) <- contr.treatment(4, base = 3)
hosp_mort_non_r3 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs_non)
summary(hosp_mort_non_r3)
##
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp +
## APACHE + non_white + copd + chf, family = binomial(link = "logit"),
## data = imgs_non)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4147 -0.6999 -0.4993 -0.2650 2.5139
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.615833 0.697808 -3.749 0.000178 ***
## quartile1 -1.041093 0.382233 -2.724 0.006455 **
## quartile2 -0.618918 0.349294 -1.772 0.076409 .
## quartile4 -0.177360 0.323987 -0.547 0.584083
## Male1 -0.114983 0.253055 -0.454 0.649556
## Age 0.007204 0.008351 0.863 0.388323
## Immsupp1 0.918398 0.259996 3.532 0.000412 ***
## APACHE 0.047948 0.016177 2.964 0.003037 **
## non_white1 -0.843357 0.373297 -2.259 0.023870 *
## copd1 -0.123341 0.346305 -0.356 0.721718
## chf1 -0.862727 0.525182 -1.643 0.100441
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 451.96 on 435 degrees of freedom
## Residual deviance: 398.98 on 425 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 420.98
##
## Number of Fisher Scoring iterations: 5
cbind(OR = exp(coef(hosp_mort_non_r3)),
exp(confint(hosp_mort_non_r3)),
p = coef(summary(hosp_mort_non_r3))[,4])
## OR 2.5 % 97.5 % p
## (Intercept) 0.07310687 0.01786178 0.2773064 0.0001777939
## quartile1 0.35306868 0.16241444 0.7333038 0.0064552941
## quartile2 0.53852667 0.26818985 1.0605029 0.0764089248
## quartile4 0.83747791 0.44198820 1.5796181 0.5840826873
## Male1 0.89138126 0.54238313 1.4658950 0.6495556183
## Age 1.00723006 0.99100117 1.0240774 0.3883234008
## Immsupp1 2.50527300 1.50706005 4.1853189 0.0004118774
## APACHE 1.04911626 1.01656587 1.0833334 0.0030366608
## non_white1 0.43026373 0.19659752 0.8615848 0.0238701170
## copd1 0.88396182 0.43715709 1.7118293 0.7217175667
## chf1 0.42200962 0.13462018 1.0960007 0.1004410074
#ref level = 4
contrasts(imgs_non$quartile) <- contr.treatment(4, base = 4)
hosp_mort_non_r4 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs_non)
summary(hosp_mort_non_r4)
##
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp +
## APACHE + non_white + copd + chf, family = binomial(link = "logit"),
## data = imgs_non)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4147 -0.6999 -0.4993 -0.2650 2.5139
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.793193 0.715721 -3.903 9.52e-05 ***
## quartile1 -0.863732 0.393273 -2.196 0.028073 *
## quartile2 -0.441558 0.358730 -1.231 0.218363
## quartile3 0.177360 0.323987 0.547 0.584083
## Male1 -0.114983 0.253055 -0.454 0.649556
## Age 0.007204 0.008351 0.863 0.388323
## Immsupp1 0.918398 0.259996 3.532 0.000412 ***
## APACHE 0.047948 0.016177 2.964 0.003037 **
## non_white1 -0.843357 0.373297 -2.259 0.023870 *
## copd1 -0.123341 0.346305 -0.356 0.721718
## chf1 -0.862727 0.525182 -1.643 0.100441
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 451.96 on 435 degrees of freedom
## Residual deviance: 398.98 on 425 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 420.98
##
## Number of Fisher Scoring iterations: 5
cbind(OR = exp(coef(hosp_mort_non_r4)),
exp(confint(hosp_mort_non_r4)),
p = coef(summary(hosp_mort_non_r4))[,4])
## OR 2.5 % 97.5 % p
## (Intercept) 0.06122539 0.01437719 0.2395434 9.515319e-05
## quartile1 0.42158566 0.19041208 0.8979589 2.807291e-02
## quartile2 0.64303388 0.31510843 1.2933171 2.183629e-01
## quartile3 1.19406135 0.63306442 2.2625038 5.840827e-01
## Male1 0.89138126 0.54238313 1.4658950 6.495556e-01
## Age 1.00723006 0.99100117 1.0240774 3.883234e-01
## Immsupp1 2.50527300 1.50706005 4.1853189 4.118774e-04
## APACHE 1.04911626 1.01656587 1.0833334 3.036661e-03
## non_white1 0.43026373 0.19659752 0.8615848 2.387012e-02
## copd1 0.88396182 0.43715709 1.7118293 7.217176e-01
## chf1 0.42200962 0.13462018 1.0960007 1.004410e-01
#reset ref level to 1
contrasts(imgs_non$quartile) <- contr.treatment(4, base = 1)
ICU Length of stay
contrasts(imgs$quartile) <-contr.treatment(4,base = 1)
#Negative binomial model
#std errors checked -> okay
ICUdays_nb <- glm.nb(ICU_days ~ Total+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=imgs)
summary(ICUdays_nb)
##
## Call:
## glm.nb(formula = ICU_days ~ Total + APACHE + Age + non_white +
## Male + Immsupp + copd + chf, data = imgs, init.theta = 3.02729853,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2276 -0.8067 -0.3650 0.2561 3.6493
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.493862 0.148292 10.074 < 2e-16 ***
## Total 0.056503 0.006971 8.106 5.25e-16 ***
## APACHE 0.022434 0.003631 6.178 6.48e-10 ***
## Age -0.009249 0.001905 -4.856 1.20e-06 ***
## non_white1 -0.242192 0.075780 -3.196 0.00139 **
## Male1 0.100367 0.059707 1.681 0.09277 .
## Immsupp1 -0.074787 0.063722 -1.174 0.24054
## copd1 -0.014389 0.084868 -0.170 0.86536
## chf1 -0.149653 0.110395 -1.356 0.17522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3.0273) family taken to be 1)
##
## Null deviance: 690.58 on 558 degrees of freedom
## Residual deviance: 533.76 on 550 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 3117.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3.027
## Std. Err.: 0.247
##
## 2 x log-likelihood: -3097.171
cbind(IRR = exp(coef(ICUdays_nb)), exp(confint(ICUdays_nb)))[2,]
## IRR 2.5 % 97.5 %
## 1.058129 1.043567 1.072946
#standard errors checked -> okay
#Base = 1
contrasts(imgs$quartile) <-contr.treatment(4,base = 1)
ICUdays_nb_q1 <- glm.nb(ICU_days ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=imgs)
summary(ICUdays_nb_q1)
##
## Call:
## glm.nb(formula = ICU_days ~ quartile + APACHE + Age + non_white +
## Male + Immsupp + copd + chf, data = imgs, init.theta = 3.025792556,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3082 -0.8335 -0.3579 0.2777 3.8017
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.516790 0.148598 10.207 < 2e-16 ***
## quartile2 0.334510 0.089146 3.752 0.000175 ***
## quartile3 0.582179 0.087980 6.617 3.66e-11 ***
## quartile4 0.654711 0.087890 7.449 9.39e-14 ***
## APACHE 0.022001 0.003651 6.026 1.68e-09 ***
## Age -0.009538 0.001923 -4.960 7.03e-07 ***
## non_white1 -0.245980 0.075787 -3.246 0.001172 **
## Male1 0.077991 0.059951 1.301 0.193288
## Immsupp1 -0.076519 0.063748 -1.200 0.230009
## copd1 -0.030197 0.085383 -0.354 0.723593
## chf1 -0.139401 0.110515 -1.261 0.207174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3.0258) family taken to be 1)
##
## Null deviance: 690.34 on 558 degrees of freedom
## Residual deviance: 533.19 on 548 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 3120.8
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3.026
## Std. Err.: 0.247
##
## 2 x log-likelihood: -3096.780
cbind(IRR = exp(coef(ICUdays_nb_q1)),
exp(confint(ICUdays_nb_q1)),
p = coef(summary(ICUdays_nb_q1))[,4])
## IRR 2.5 % 97.5 % p
## (Intercept) 4.5575713 3.3769015 6.1566319 1.837967e-24
## quartile2 1.3972550 1.1729172 1.6648426 1.751569e-04
## quartile3 1.7899342 1.5048166 2.1296709 3.660726e-11
## quartile4 1.9245869 1.6188460 2.2888263 9.389953e-14
## APACHE 1.0222452 1.0146871 1.0298937 1.676954e-09
## Age 0.9905074 0.9866450 0.9943669 7.032483e-07
## non_white1 0.7819378 0.6742871 0.9077580 1.171743e-03
## Male1 1.0811128 0.9606191 1.2165457 1.932881e-01
## Immsupp1 0.9263357 0.8178175 1.0494715 2.300089e-01
## copd1 0.9702548 0.8213731 1.1482057 7.235930e-01
## chf1 0.8698788 0.7014168 1.0830270 2.071737e-01
#Base = 2
contrasts(imgs$quartile) <-contr.treatment(4,base = 2)
ICUdays_nb_q2 <- glm.nb(ICU_days ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=imgs)
summary(ICUdays_nb_q2)
##
## Call:
## glm.nb(formula = ICU_days ~ quartile + APACHE + Age + non_white +
## Male + Immsupp + copd + chf, data = imgs, init.theta = 3.025792556,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3082 -0.8335 -0.3579 0.2777 3.8017
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.851299 0.161818 11.441 < 2e-16 ***
## quartile1 -0.334510 0.089146 -3.752 0.000175 ***
## quartile3 0.247669 0.083791 2.956 0.003119 **
## quartile4 0.320202 0.083407 3.839 0.000124 ***
## APACHE 0.022001 0.003651 6.026 1.68e-09 ***
## Age -0.009538 0.001923 -4.960 7.03e-07 ***
## non_white1 -0.245980 0.075787 -3.246 0.001172 **
## Male1 0.077991 0.059951 1.301 0.193288
## Immsupp1 -0.076519 0.063748 -1.200 0.230009
## copd1 -0.030197 0.085383 -0.354 0.723593
## chf1 -0.139401 0.110515 -1.261 0.207174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3.0258) family taken to be 1)
##
## Null deviance: 690.34 on 558 degrees of freedom
## Residual deviance: 533.19 on 548 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 3120.8
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3.026
## Std. Err.: 0.247
##
## 2 x log-likelihood: -3096.780
cbind(IRR = exp(coef(ICUdays_nb_q2)),
exp(confint(ICUdays_nb_q2)),
p = coef(summary(ICUdays_nb_q2))[,4])
## IRR 2.5 % 97.5 % p
## (Intercept) 6.3680892 4.5987356 8.8315486 2.620331e-30
## quartile1 0.7156890 0.6006574 0.8525751 1.751569e-04
## quartile3 1.2810362 1.0871930 1.5095506 3.118715e-03
## quartile4 1.3774056 1.1684231 1.6239800 1.235251e-04
## APACHE 1.0222452 1.0146871 1.0298937 1.676954e-09
## Age 0.9905074 0.9866450 0.9943669 7.032483e-07
## non_white1 0.7819378 0.6742871 0.9077580 1.171743e-03
## Male1 1.0811128 0.9606191 1.2165457 1.932881e-01
## Immsupp1 0.9263357 0.8178175 1.0494715 2.300089e-01
## copd1 0.9702548 0.8213731 1.1482057 7.235930e-01
## chf1 0.8698788 0.7014168 1.0830270 2.071737e-01
#Base = 3
contrasts(imgs$quartile) <-contr.treatment(4,base = 3)
ICUdays_nb_q3 <- glm.nb(ICU_days ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=imgs)
summary(ICUdays_nb_q3)
##
## Call:
## glm.nb(formula = ICU_days ~ quartile + APACHE + Age + non_white +
## Male + Immsupp + copd + chf, data = imgs, init.theta = 3.025792556,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3082 -0.8335 -0.3579 0.2777 3.8017
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.098969 0.157668 13.313 < 2e-16 ***
## quartile1 -0.582179 0.087980 -6.617 3.66e-11 ***
## quartile2 -0.247669 0.083791 -2.956 0.00312 **
## quartile4 0.072532 0.081171 0.894 0.37155
## APACHE 0.022001 0.003651 6.026 1.68e-09 ***
## Age -0.009538 0.001923 -4.960 7.03e-07 ***
## non_white1 -0.245980 0.075787 -3.246 0.00117 **
## Male1 0.077991 0.059951 1.301 0.19329
## Immsupp1 -0.076519 0.063748 -1.200 0.23001
## copd1 -0.030197 0.085383 -0.354 0.72359
## chf1 -0.139401 0.110515 -1.261 0.20717
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3.0258) family taken to be 1)
##
## Null deviance: 690.34 on 558 degrees of freedom
## Residual deviance: 533.19 on 548 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 3120.8
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3.026
## Std. Err.: 0.247
##
## 2 x log-likelihood: -3096.780
cbind(IRR = exp(coef(ICUdays_nb_q3)),
exp(confint(ICUdays_nb_q3)),
p = coef(summary(ICUdays_nb_q3))[,4])
## IRR 2.5 % 97.5 % p
## (Intercept) 8.1577527 5.9060823 11.2873406 1.955984e-40
## quartile1 0.5586798 0.4695561 0.6645328 3.660726e-11
## quartile2 0.7806181 0.6624488 0.9197999 3.118715e-03
## quartile4 1.0752277 0.9167725 1.2611534 3.715475e-01
## APACHE 1.0222452 1.0146871 1.0298937 1.676954e-09
## Age 0.9905074 0.9866450 0.9943669 7.032483e-07
## non_white1 0.7819378 0.6742871 0.9077580 1.171743e-03
## Male1 1.0811128 0.9606191 1.2165457 1.932881e-01
## Immsupp1 0.9263357 0.8178175 1.0494715 2.300089e-01
## copd1 0.9702548 0.8213731 1.1482057 7.235930e-01
## chf1 0.8698788 0.7014168 1.0830270 2.071737e-01
#Base = 4
contrasts(imgs$quartile) <-contr.treatment(4,base = 4)
ICUdays_nb_q4 <- glm.nb(ICU_days ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=imgs)
summary(ICUdays_nb_q4)
##
## Call:
## glm.nb(formula = ICU_days ~ quartile + APACHE + Age + non_white +
## Male + Immsupp + copd + chf, data = imgs, init.theta = 3.025792556,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3082 -0.8335 -0.3579 0.2777 3.8017
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.171501 0.159360 13.626 < 2e-16 ***
## quartile1 -0.654711 0.087890 -7.449 9.39e-14 ***
## quartile2 -0.320202 0.083407 -3.839 0.000124 ***
## quartile3 -0.072532 0.081171 -0.894 0.371548
## APACHE 0.022001 0.003651 6.026 1.68e-09 ***
## Age -0.009538 0.001923 -4.960 7.03e-07 ***
## non_white1 -0.245980 0.075787 -3.246 0.001172 **
## Male1 0.077991 0.059951 1.301 0.193288
## Immsupp1 -0.076519 0.063748 -1.200 0.230009
## copd1 -0.030197 0.085383 -0.354 0.723593
## chf1 -0.139401 0.110515 -1.261 0.207174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3.0258) family taken to be 1)
##
## Null deviance: 690.34 on 558 degrees of freedom
## Residual deviance: 533.19 on 548 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 3120.8
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3.026
## Std. Err.: 0.247
##
## 2 x log-likelihood: -3096.780
cbind(IRR = exp(coef(ICUdays_nb_q4)),
exp(confint(ICUdays_nb_q4)),
p = coef(summary(ICUdays_nb_q4))[,4])
## IRR 2.5 % 97.5 % p
## (Intercept) 8.7714418 6.3617496 12.1155580 2.790730e-42
## quartile1 0.5195920 0.4369051 0.6177240 9.389953e-14
## quartile2 0.7260026 0.6157711 0.8558543 1.235251e-04
## quartile3 0.9300355 0.7929249 1.0907832 3.715475e-01
## APACHE 1.0222452 1.0146871 1.0298937 1.676954e-09
## Age 0.9905074 0.9866450 0.9943669 7.032483e-07
## non_white1 0.7819378 0.6742871 0.9077580 1.171743e-03
## Male1 1.0811128 0.9606191 1.2165457 1.932881e-01
## Immsupp1 0.9263357 0.8178175 1.0494715 2.300089e-01
## copd1 0.9702548 0.8213731 1.1482057 7.235930e-01
## chf1 0.8698788 0.7014168 1.0830270 2.071737e-01
#Reset base level
contrasts(imgs$quartile) <-contr.treatment(4,base = 1)
Duration of mechanical ventilation
#make the subset of data: only intubated patients (n = 286)
intub <- imgs %>% filter(days_MV != 0)
#standard errors checked ->okay
MVdays_nb2 <- glm.nb(days_MV ~ Total+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=intub)
summary(MVdays_nb2)
##
## Call:
## glm.nb(formula = days_MV ~ Total + APACHE + Age + non_white +
## Male + Immsupp + copd + chf, data = intub, init.theta = 1.805548498,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7943 -0.9164 -0.4188 0.2213 4.3968
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.139983 0.260729 8.208 2.26e-16 ***
## Total 0.044079 0.012079 3.649 0.000263 ***
## APACHE 0.003229 0.006372 0.507 0.612389
## Age -0.008990 0.003171 -2.835 0.004576 **
## non_white1 0.016122 0.137941 0.117 0.906958
## Male1 0.093835 0.099917 0.939 0.347663
## Immsupp1 -0.021805 0.104509 -0.209 0.834729
## copd1 -0.238111 0.145005 -1.642 0.100573
## chf1 0.041066 0.196049 0.209 0.834081
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.8055) family taken to be 1)
##
## Null deviance: 315.11 on 285 degrees of freedom
## Residual deviance: 290.28 on 277 degrees of freedom
## AIC: 1776.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.806
## Std. Err.: 0.174
##
## 2 x log-likelihood: -1756.851
cbind(IRR = exp(coef(MVdays_nb2)), exp(confint(MVdays_nb2)))[2,]
## IRR 2.5 % 97.5 %
## 1.045064 1.020284 1.070437
#Standard errors checked -> okay
MVdays_nbq <- glm.nb(days_MV ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=intub)
summary(MVdays_nbq)
##
## Call:
## glm.nb(formula = days_MV ~ quartile + APACHE + Age + non_white +
## Male + Immsupp + copd + chf, data = intub, init.theta = 1.828470791,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8183 -0.9216 -0.4349 0.2329 4.2200
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.212944 0.262174 8.441 < 2e-16 ***
## quartile2 0.020217 0.167878 0.120 0.90415
## quartile3 0.434000 0.161212 2.692 0.00710 **
## quartile4 0.449699 0.157266 2.859 0.00424 **
## APACHE 0.001820 0.006392 0.285 0.77586
## Age -0.008330 0.003165 -2.632 0.00849 **
## non_white1 0.015970 0.137671 0.116 0.90765
## Male1 0.096665 0.100921 0.958 0.33815
## Immsupp1 -0.027713 0.104094 -0.266 0.79006
## copd1 -0.275100 0.145398 -1.892 0.05848 .
## chf1 0.029604 0.195470 0.151 0.87962
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.8285) family taken to be 1)
##
## Null deviance: 318.28 on 285 degrees of freedom
## Residual deviance: 289.70 on 275 degrees of freedom
## AIC: 1777.4
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.828
## Std. Err.: 0.176
##
## 2 x log-likelihood: -1753.373
#cbind(IRR = exp(coef(MVdays_nbq)), exp(confint(MVdays_nbq)))[2:4,]
#Base = 1
contrasts(intub$quartile) <-contr.treatment(4,base = 1)
MVdays_nb2_q1 <- glm.nb(days_MV ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=intub)
summary(MVdays_nb2_q1)
##
## Call:
## glm.nb(formula = days_MV ~ quartile + APACHE + Age + non_white +
## Male + Immsupp + copd + chf, data = intub, init.theta = 1.828470791,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8183 -0.9216 -0.4349 0.2329 4.2200
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.212944 0.262174 8.441 < 2e-16 ***
## quartile2 0.020217 0.167878 0.120 0.90415
## quartile3 0.434000 0.161212 2.692 0.00710 **
## quartile4 0.449699 0.157266 2.859 0.00424 **
## APACHE 0.001820 0.006392 0.285 0.77586
## Age -0.008330 0.003165 -2.632 0.00849 **
## non_white1 0.015970 0.137671 0.116 0.90765
## Male1 0.096665 0.100921 0.958 0.33815
## Immsupp1 -0.027713 0.104094 -0.266 0.79006
## copd1 -0.275100 0.145398 -1.892 0.05848 .
## chf1 0.029604 0.195470 0.151 0.87962
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.8285) family taken to be 1)
##
## Null deviance: 318.28 on 285 degrees of freedom
## Residual deviance: 289.70 on 275 degrees of freedom
## AIC: 1777.4
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.828
## Std. Err.: 0.176
##
## 2 x log-likelihood: -1753.373
cbind(IRR = exp(coef(MVdays_nb2_q1)),
exp(confint(MVdays_nb2_q1)),
p = coef(summary(MVdays_nb2_q1))[,4])
## IRR 2.5 % 97.5 % p
## (Intercept) 9.1425923 5.3009692 15.8962705 3.153570e-17
## quartile2 1.0204223 0.7272564 1.4256791 9.041472e-01
## quartile3 1.5434189 1.1091160 2.1368015 7.100218e-03
## quartile4 1.5678399 1.1401674 2.1435081 4.243256e-03
## APACHE 1.0018216 0.9886235 1.0151979 7.758557e-01
## Age 0.9917047 0.9851211 0.9982833 8.489914e-03
## non_white1 1.0160986 0.7766884 1.3417464 9.076497e-01
## Male1 1.1014916 0.8993489 1.3480047 3.381473e-01
## Immsupp1 0.9726676 0.7968709 1.1890068 7.900622e-01
## copd1 0.7594960 0.5735750 1.0140200 5.848437e-02
## chf1 1.0300465 0.7087438 1.5327841 8.796204e-01
#Base = 2
contrasts(intub$quartile) <-contr.treatment(4,base = 2)
MVdays_nb2_q2 <- glm.nb(days_MV ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=intub)
summary(MVdays_nb2_q2)
##
## Call:
## glm.nb(formula = days_MV ~ quartile + APACHE + Age + non_white +
## Male + Immsupp + copd + chf, data = intub, init.theta = 1.828470791,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8183 -0.9216 -0.4349 0.2329 4.2200
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.233161 0.274568 8.133 4.18e-16 ***
## quartile1 -0.020217 0.167878 -0.120 0.90415
## quartile3 0.413783 0.137860 3.001 0.00269 **
## quartile4 0.429482 0.134143 3.202 0.00137 **
## APACHE 0.001820 0.006392 0.285 0.77586
## Age -0.008330 0.003165 -2.632 0.00849 **
## non_white1 0.015970 0.137671 0.116 0.90765
## Male1 0.096665 0.100921 0.958 0.33815
## Immsupp1 -0.027713 0.104094 -0.266 0.79006
## copd1 -0.275100 0.145398 -1.892 0.05848 .
## chf1 0.029604 0.195470 0.151 0.87962
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.8285) family taken to be 1)
##
## Null deviance: 318.28 on 285 degrees of freedom
## Residual deviance: 289.70 on 275 degrees of freedom
## AIC: 1777.4
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.828
## Std. Err.: 0.176
##
## 2 x log-likelihood: -1753.373
cbind(IRR = exp(coef(MVdays_nb2_q2)),
exp(confint(MVdays_nb2_q2)),
p = coef(summary(MVdays_nb2_q2))[,4])
## IRR 2.5 % 97.5 % p
## (Intercept) 9.3293051 5.3385349 16.4324192 4.175382e-16
## quartile1 0.9799864 0.7014201 1.3750308 9.041472e-01
## quartile3 1.5125295 1.1496822 1.9884426 2.686837e-03
## quartile4 1.5364618 1.1785024 2.0001430 1.366298e-03
## APACHE 1.0018216 0.9886235 1.0151979 7.758557e-01
## Age 0.9917047 0.9851211 0.9982833 8.489914e-03
## non_white1 1.0160986 0.7766884 1.3417464 9.076497e-01
## Male1 1.1014916 0.8993489 1.3480047 3.381473e-01
## Immsupp1 0.9726676 0.7968709 1.1890068 7.900622e-01
## copd1 0.7594960 0.5735750 1.0140200 5.848437e-02
## chf1 1.0300465 0.7087438 1.5327841 8.796204e-01
#Base = 3
contrasts(intub$quartile) <-contr.treatment(4,base = 3)
MVdays_nb2_q3 <- glm.nb(days_MV ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=intub)
summary(MVdays_nb2_q3)
##
## Call:
## glm.nb(formula = days_MV ~ quartile + APACHE + Age + non_white +
## Male + Immsupp + copd + chf, data = intub, init.theta = 1.828470791,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8183 -0.9216 -0.4349 0.2329 4.2200
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.646944 0.269969 9.805 < 2e-16 ***
## quartile1 -0.434000 0.161212 -2.692 0.00710 **
## quartile2 -0.413783 0.137860 -3.001 0.00269 **
## quartile4 0.015699 0.123956 0.127 0.89922
## APACHE 0.001820 0.006392 0.285 0.77586
## Age -0.008330 0.003165 -2.632 0.00849 **
## non_white1 0.015970 0.137671 0.116 0.90765
## Male1 0.096665 0.100921 0.958 0.33815
## Immsupp1 -0.027713 0.104094 -0.266 0.79006
## copd1 -0.275100 0.145398 -1.892 0.05848 .
## chf1 0.029604 0.195470 0.151 0.87962
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.8285) family taken to be 1)
##
## Null deviance: 318.28 on 285 degrees of freedom
## Residual deviance: 289.70 on 275 degrees of freedom
## AIC: 1777.4
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.828
## Std. Err.: 0.176
##
## 2 x log-likelihood: -1753.373
cbind(IRR = exp(coef(MVdays_nb2_q3)),
exp(confint(MVdays_nb2_q3)),
p = coef(summary(MVdays_nb2_q3))[,4])
## IRR 2.5 % 97.5 % p
## (Intercept) 14.1108494 7.9937450 25.0773939 1.075494e-22
## quartile1 0.6479123 0.4679892 0.9016189 7.100218e-03
## quartile2 0.6611441 0.5029062 0.8698056 2.686837e-03
## quartile4 1.0158227 0.7939350 1.2988435 8.992187e-01
## APACHE 1.0018216 0.9886235 1.0151979 7.758557e-01
## Age 0.9917047 0.9851211 0.9982833 8.489914e-03
## non_white1 1.0160986 0.7766884 1.3417464 9.076497e-01
## Male1 1.1014916 0.8993489 1.3480047 3.381473e-01
## Immsupp1 0.9726676 0.7968709 1.1890068 7.900622e-01
## copd1 0.7594960 0.5735750 1.0140200 5.848437e-02
## chf1 1.0300465 0.7087438 1.5327841 8.796204e-01
#Base = 4
contrasts(intub$quartile) <-contr.treatment(4,base = 4)
MVdays_nb2_q4 <- glm.nb(days_MV ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=intub)
summary(MVdays_nb2_q4)
##
## Call:
## glm.nb(formula = days_MV ~ quartile + APACHE + Age + non_white +
## Male + Immsupp + copd + chf, data = intub, init.theta = 1.828470791,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8183 -0.9216 -0.4349 0.2329 4.2200
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.662643 0.263364 10.110 < 2e-16 ***
## quartile1 -0.449699 0.157266 -2.859 0.00424 **
## quartile2 -0.429482 0.134143 -3.202 0.00137 **
## quartile3 -0.015699 0.123956 -0.127 0.89922
## APACHE 0.001820 0.006392 0.285 0.77586
## Age -0.008330 0.003165 -2.632 0.00849 **
## non_white1 0.015970 0.137671 0.116 0.90765
## Male1 0.096665 0.100921 0.958 0.33815
## Immsupp1 -0.027713 0.104094 -0.266 0.79006
## copd1 -0.275100 0.145398 -1.892 0.05848 .
## chf1 0.029604 0.195470 0.151 0.87962
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.8285) family taken to be 1)
##
## Null deviance: 318.28 on 285 degrees of freedom
## Residual deviance: 289.70 on 275 degrees of freedom
## AIC: 1777.4
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.828
## Std. Err.: 0.176
##
## 2 x log-likelihood: -1753.373
cbind(IRR = exp(coef(MVdays_nb2_q4)),
exp(confint(MVdays_nb2_q4)),
p = coef(summary(MVdays_nb2_q4))[,4])
## IRR 2.5 % 97.5 % p
## (Intercept) 14.3341208 8.3166527 24.8848017 4.981497e-24
## quartile1 0.6378202 0.4665249 0.8770642 4.243256e-03
## quartile2 0.6508460 0.4999643 0.8485345 1.366298e-03
## quartile3 0.9844238 0.7699157 1.2595489 8.992187e-01
## APACHE 1.0018216 0.9886235 1.0151979 7.758557e-01
## Age 0.9917047 0.9851211 0.9982833 8.489914e-03
## non_white1 1.0160986 0.7766884 1.3417464 9.076497e-01
## Male1 1.1014916 0.8993489 1.3480047 3.381473e-01
## Immsupp1 0.9726676 0.7968709 1.1890068 7.900622e-01
## copd1 0.7594960 0.5735750 1.0140200 5.848437e-02
## chf1 1.0300465 0.7087438 1.5327841 8.796204e-01
#Reset base level
contrasts(intub$quartile) <-contr.treatment(4,base = 1)
#ICU days (mean x IRR)
incr_days = (1.058129 * mean(imgs$ICU_days)) - mean(imgs$ICU_days)
incr_days
## [1] 0.4002597
median(imgs$ICU_days)
## [1] 5
#Days MV
more_days = (1.045064 * mean(intub$days_MV)) - mean(intub$days_MV)
more_days
## [1] 0.3691782
median(intub$days_MV)
## [1] 5
library(reshape2)
library(RColorBrewer)
#Eliminate non-intubated from the MV data and transform for plotting
meltICU <- imgs %>% dplyr::select(Subject_ID, ICU_days, quartile) %>% melt(id.vars="quartile", measure.vars="ICU_days")
meltMV <- imgs %>% filter(days_MV != 0) %>% dplyr::select(Subject_ID, days_MV, quartile) %>% melt(id.vars="quartile", measure.vars="days_MV")
meltplot <- rbind(meltICU, meltMV)
ICUdata <- imgs %>% dplyr::select(Subject_ID, ICU_days, quartile) %>% rename(value = ICU_days) %>% mutate(variable = "ICU_days")
MVdata <- intub %>% dplyr::select(Subject_ID, days_MV, quartile) %>% rename(value = days_MV) %>% mutate(variable = "days_MV")
meltplot <- rbind(ICUdata, MVdata)
#Boxplot square root scale
ggplot(meltplot) + geom_boxplot(aes(quartile, sqrt(value), color=variable)) +
scale_x_discrete(name="Chest Radiograph Score Quartile") +
scale_y_continuous(name="(Square root of) Days") +
theme_minimal() + theme(plot.title = element_text(face = "bold")) +
theme(panel.grid.minor = element_blank()) +
theme(legend.title = element_blank()) +
theme(legend.position = "bottom") +
scale_color_discrete(labels=c("Length of stay", "Ventilator days"))

#Boxplot log scale
ggplot(meltplot) + geom_boxplot(aes(quartile, log(value), color=variable)) +
scale_x_discrete(name="Chest Radiograph Score Quartile") +
scale_y_continuous(name="(Log of) Days", limits=c(0,5)) +
theme_minimal() + theme(plot.title = element_text(face = "bold")) +
theme(panel.grid.minor = element_blank()) +
theme(legend.title = element_blank()) +
theme(legend.position = "bottom") +
scale_color_discrete(labels=c("Length of stay", "Ventilator days"))

box_comparisons <- list(c("1","2"), c("2","3"), c("3","4"), c("1","3"), c("1","4"), c("2","4"))
symnum.args <- c("***"=0.05,ns=1)
var_names <- c('ICU_days' = "Length of ICU stay", 'days_MV' = "Ventilator days")
my_boxplot <- ggplot(meltplot, aes(quartile, log(value))) + geom_boxplot() +
facet_grid(.~ variable, labeller=as_labeller(var_names)) +
my_stat_compare_means(label = "p.signif", method = "t.test", comparisons = box_comparisons, symnum.args = symnum.args) +
scale_x_discrete(name="Chest Radiograph Score Quartile") +
scale_y_continuous(name="(Log) Days", limits=c(-0.25,7.5)) +
theme_bw() +
annotate("text", x=2.5, y=-0.2, label="p for trend < .001")
my_boxplot

wilcox_data <- meltplot %>% filter(variable == "days_MV") %>% filter(quartile == 2 | quartile == 3)
t.test(log(wilcox_data$value) ~ wilcox_data$quartile)
##
## Welch Two Sample t-test
##
## data: log(wilcox_data$value) by wilcox_data$quartile
## t = -2.0012, df = 148.09, p-value = 0.0472
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.56153943 -0.00353892
## sample estimates:
## mean in group 2 mean in group 3
## 1.512006 1.794545
# jonkheere Terpstra test
ICUdata$quartile <- ordered(ICUdata$quartile, levels=c("1", "2", "3", "4"))
jonckheere.test(ICUdata$value, ICUdata$quartile, alternative="two.sided")
## Warning in jonckheere.test(ICUdata$value, ICUdata$quartile, alternative = "two.sided"): Sample size > 100 or data with ties
## p-value based on normal approximation. Specify nperm for permutation p-value
##
## Jonckheere-Terpstra test
##
## data:
## JT = 74022, p-value = 1.16e-12
## alternative hypothesis: two.sided
MVdata$quartile <- ordered(MVdata$quartile, levels=c("1", "2", "3", "4"))
jonckheere.test(MVdata$value, MVdata$quartile, alternative="two.sided")
## Warning in jonckheere.test(MVdata$value, MVdata$quartile, alternative = "two.sided"): Sample size > 100 or data with ties
## p-value based on normal approximation. Specify nperm for permutation p-value
##
## Jonckheere-Terpstra test
##
## data:
## JT = 18238, p-value = 3.097e-05
## alternative hypothesis: two.sided
Miscellaneous
shapiro.test(imgs$survival)
##
## Shapiro-Wilk normality test
##
## data: imgs$survival
## W = 0.79791, p-value < 2.2e-16
imgs %>% summarize(mean = mean(survival), median = median(survival), min = min(survival), max = max(survival), count = length(survival))
## # A tibble: 1 x 5
## mean median min max count
## <dbl> <dbl> <dbl> <dbl> <int>
## 1 932. 364. 0 3664 560
imgs %>% mutate(surv_quart = ntile(survival, 4)) %>% group_by(surv_quart) %>% summarize(mean = mean(survival), median = median(survival), min = min(survival), max = max(survival), count = length(survival))
## # A tibble: 4 x 6
## surv_quart mean median min max count
## <int> <dbl> <dbl> <dbl> <dbl> <int>
## 1 1 8.89 8 0 22 140
## 2 2 127. 104 22 362 140
## 3 3 951. 912. 365 1788 140
## 4 4 2643. 2388 1796 3664 140
imgs %>% group_by(Final_Dx) %>% summarize(count = length(Total))
## # A tibble: 4 x 2
## Final_Dx count
## <chr> <int>
## 1 ARDS 26
## 2 Sepsis 292
## 3 Sepsis/ARDS 97
## 4 SIRS 145
imgs %>% ggplot(aes(Total, fill = quartile)) + geom_histogram(bins = 15)

imgs %>% summarize( median = median(Total), mean = mean(Total))
## # A tibble: 1 x 2
## median mean
## <dbl> <dbl>
## 1 7 6.61
imgs %>% group_by(quartile) %>% summarize(min = min(Total), max = max(Total))
## # A tibble: 4 x 3
## quartile min max
## <fct> <dbl> <dbl>
## 1 1 0 2
## 2 2 2 7
## 3 3 7 11
## 4 4 11 16
imgs %>% summarize(med_age = median(Age), med_apa = median(APACHE))
## # A tibble: 1 x 2
## med_age med_apa
## <dbl> <dbl>
## 1 60 25
imgs %>% group_by(Male) %>% summarize(count = length(Total))
## # A tibble: 2 x 2
## Male count
## <fct> <int>
## 1 0 260
## 2 1 300
imgs %>% group_by(non_white) %>% summarize(count = length(Total))
## # A tibble: 2 x 2
## non_white count
## <fct> <int>
## 1 0 437
## 2 1 123
imgs %>% group_by(Immsupp) %>% summarize(count = length(Total))
## # A tibble: 2 x 2
## Immsupp count
## <fct> <int>
## 1 0 355
## 2 1 205
imgs %>% mutate(age_qt = ntile(Age,4)) %>% group_by(age_qt) %>% summarize(lower = min(Age), upper = max(Age))
## # A tibble: 4 x 3
## age_qt lower upper
## <int> <dbl> <dbl>
## 1 1 19 47
## 2 2 47 60
## 3 3 60 70
## 4 4 70 94
imgs %>% mutate(apa_qt = ntile(APACHE,4)) %>% group_by(apa_qt) %>% summarize(lower = min(APACHE), upper = max(APACHE))
## # A tibble: 4 x 3
## apa_qt lower upper
## <int> <dbl> <dbl>
## 1 1 3 20
## 2 2 20 25
## 3 3 25 31
## 4 4 31 54
imgs %>% filter(!is.na(date_death)) %>% summarize(count = length(Total)) #total death events in cohort (Cox event rate)
## # A tibble: 1 x 1
## count
## <int>
## 1 369
imgs %>% filter(DC_coded == 0 | DC_coded == 1) %>% summarize(count = length(Total)) #in hospital deaths - entire cohort
## # A tibble: 1 x 1
## count
## <int>
## 1 159
imgs_non %>% filter(DC_coded == 0 | DC_coded == 1) %>% summarize(count = length(Total)) # in hospital deaths for ARDS only
## # A tibble: 1 x 1
## count
## <int>
## 1 93
key <- read_excel("~/Documents/ROCIanalysis/ROCI_101918.xlsx", sheet = "Key", range = cell_cols(1:2))
subjects <- key %>% filter(Subject_ID %in% imgs$Subject_ID)
subjects <- imgs %>% dplyr::select(Subject_ID, cxr_date, Hosp_dc) %>% right_join(subjects)
write.csv(subjects, file="subjects.csv", row.names = FALSE)
#CXR scores are non-normally distributed and ordinal, will use a spearman correlation.
imgs %>% dplyr::select(Total, APACHE) %>% cor(., use = "pairwise.complete.obs", method="spearman")
## Total APACHE
## Total 1.0000000 0.1764567
## APACHE 0.1764567 1.0000000
cor.test(imgs$Total, imgs$APACHE, method="spearman", alternative = "t")
## Warning in cor.test.default(imgs$Total, imgs$APACHE, method = "spearman", :
## Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: imgs$Total and imgs$APACHE
## S = 24104000, p-value = 2.676e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.1764567
imgs %>% ggplot(aes(Subject_ID, Total, color=hosp_death)) + geom_jitter() + scale_color_manual(values=c("black", "red"))

imgs %>% ggplot(aes(quartile, Total, color=hosp_death)) + geom_jitter() + scale_color_manual(values=c("black", "red"))

imgs %>% ggplot(aes(Total, fill = hosp_death)) + geom_histogram() + scale_fill_manual(values=c("black", "red"))
